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ABSTRACT 



The partial differential equation for heat diffusion is 
numerically Integrated by the Runga-Kutta method. Solutions 
are obtained for the diurnal temperature variation with a 
bounded coefficient of eddy diffusivity which varies with the 
lapse rate and with height. The surface wave is represented 
by the sum of a diurnal and a semidiurnal harmonic wave. The 
results may be interpreted to apply over a fairly broad range 
of diffusivity with height. With appropriate choices of the 
various parameters, reasonably good agreement is obtained be* 
tween theoretical and observational values of amplitude and phase 
lag as functions of height and time. 
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Introduction. 



Many efforts have been made by meteorologists to solve 
the heat diffusion problem in the atmosphere, in particular the 
diurnal temperature wave has received much attention. In every 
attempt the concept of eddy diffusivity has been retained, the 
various functional forms of eddy diffusivity give different 
results. Most analyses are based upon the solution of the Taylor 
heat diffusion equation, which states that the turbulent flux of 
heat is proportional to the gradient of potential temperature, Namely 




Here t represents the time, H. the height, and the potential 

temperature deviation from a constant (mean) value. 

The simplest approach is based upon the assumption of a 
constant diffusion coefficient K. According to this solution the 
amplitude of the diurnal temperature wave decreases exponentially 
and phase increases linearly with height. For a variable K 
several solutions have been presented. Sutton (1) in his compre- 
hensive survey of daily temperature variation has considered K - Kj,z 
where ra is the stability parameter having magnitude O yr\. \ . 

His final solution is in terras of Bessel functions. It may be 

inferred from his results that at sufficiently great heights, the 

1 1 

diurnal variation becomes periodic with a phase lag which is 
proportional to z , so that the exponent of z is now as 

low as 0.5 (for m=l) . However, agreement with observation is not as 
good as might be expected in view of the additional complexity of 
the solution. 
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Best and Kohler (2) have made a similar approach to the 
problem. Haurwitz (3) has considered K as: ko Ko + K^z , where 
Ko and R. are constants. By proper selection of Ko and K. 

* * t 

Haurwitz obtains good agreement with data taken at Leafield, 
England, up to 10 meters. In his results, the amplitude falls 
off more rapidly at low levels than in K Constant Theory. How- 
ever, it is greater at high levels where K is larger. The phase 
lag is larger near the ground where K is small but becomes smaller 
at higher levels where K is larger. Poppendiek (4) has obtained 
a solution, for K varying sinusoidally with time and linearly 
with height-. The solution is very complicated, consisting of an 
infinite trigonometric series having coefficients determined by 
difficult integrations of functions of the real and imaginary 
parts of Hankel functions of complex argument. No numerical values 
are given, and the results are not suitable for practical purposes. 
Even for the case of K Independent of time but varying as a power 
of height, the solution is in terms of Bessel functions of order 
depending on this power. Staley (5) has investigated the' problem 
with K increasing with height but bounded. In several respects 
the results show better agreement with observations than previous 
solutions for the coefficient unbounded. A numerical solution 
with a diffusion coefficient which varies periodically with time 
and exponentially with height has been obtained by Haltiner (7). 

In general the functional forms of the diffusion coefficient have 
hot included such parameters as the roughness coefficient, gradient 
wind, lapse rate etc., which are generally believed to be closely 
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related to the diffusion process. It is the purpose of this investiga- 
tion to attempt to determine a diffusion coefficient in which the time 
variation is related to the stability. Briefly, the present analysis 
consists of the selection of proper functional form for the coefficient 
of eddy diffusivity from observed data, followed by numerical solution 
of the Taylor heat diffusion equation for this value of K. 
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2. Heat Dlffusivlty Coefficient. 

From physical considerations, there is no reason to expect 
an ever increasing K. Hence a bounded function of K appears to 
be a logical choice. Also, many studies on the variation of K 
with various meteorological parameters give clues to the form- 
ation of a suitable function. Sutton (1) has remarked that " the 
relatively small change of amplitude of the diurnal temperature 

a 

wave above 10 meters in summer compared with that in winter 
indicates that a more effective mechanism for the upward transfer 
of heat is at work in the summer than in the winter". Also 
Haltiner (7) suggests that K should have diurnal variation. 

It is possible that both the seasonal and diurnal variation in 
diffusion may reasonably be accounted for by making K a suitable 
function of lapse rate. Moreover from the observations made at the 
micrometeorological tower at the University of Washington, Fig. (1), 
which shows K as a function of time at three heights, it is evident 
ttiat K at a given level increases rapidly by one or more orders of 
magnitude as the lapse rate changes from sub-adiabatic to super- 
adiabatic after - sunrise. In the evening when the lapse rate changes 
from super-adiabatic to slightly sub-adiabatic K drops rapidly by one 
or more order of magnitude at any given level. Further from Fig. (2), 
which (hows K as a function of lapse rate, it is evident that at lower 
levels the variation in K is less than at higher levels. On further 
investigation of the data collected at the University of Texas micro- 
meteorological tower (see Table 1) it has been found that the variation 
in K at the lower levels is an order of magnitude (factor of ten) while 
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V 



Fig. I 




Coefficient of eddy conductivity as a function of time at 
different heights; a: 0.75 m at University of Washington micro* 
meteorological mast and tower, represents mean over 23 days of 
Aug. 1950, chosen for having most nearly sinusoidal temperature 
variations; b and c : for 15 and 35 m, respectively ( after Jehn 
and Gerherdt, 1950) ,3*4 Aug 1948, at Manor Texas. Dashed 
lines represent unreliability. All values obtained by methods 
based on energy continuity at earths surface. 
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at higher levels K varies by two orders of magnitude. Therefore 
any function which represents K should show the above character- 
istics. As a preliminary study the following functions have been 
investigated: 



K = 


a -[' -0- t + a 5 


•***! 


(2) 




f - a 3 






K r 


a , - a i e ' 


M J 


( 3 ) 


K - 


a, (| 

IZ 1 




W 


K - 


a.,|j - do 


JH^'J 


( 5 ) 


represents the height. 


0 is the potential 


temperature 



deviation from mean value and a^» a^ » a^ , a , a^ are contents. 

In equation (2) the second factor ( ' - a + e * a ‘ 5 ' 1 ) is obviously 
bounded; and the other factor (j + changes K as lapse 

rate changes. Furthermore, the variable z is also involved in this 
factor. This is necessary, because at great heights the change in 
lapse rate becomes very small and the factor z is needed to amplify 
the effect of lapse rate at these elevations. When z » 0, bg. is 
maximum, then K reduces to K ■ ) ( i — a.^ ) for z ■ oO , ^5. - o & 

KrCL|In quality this equation possesses all characteristics mentioned 
above but the quantitative agreement with the observed data is not 
sufficiently good. In the second representation (3) K is bounded 
and it increases with height to a certain level which may be deter- 
mined by the proper selection of constants. There is an appreciable 
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difference in day-time and night-time, but the night-time values are lower 

at high elevations than at lower levels. This function is very sensitive 

0 

to changes of and in magnitude of z. The square of the quantity 

was intended to avoid negative values of K. 

Form (4) is somewhat simpler than equation (3), In this equation also, 

K is bounded and has much the same characteristics as equation (3), but the 
variation in the value of K when lapse rate changes sign is less than the 
former. Diffusion coefficient (5) gives a reasonably good fit with the 
observed data although the night-time values are higher than observed. 

The difference between night and day-time values is correctly maintained. 

For the present analysis, equation (5) has been selected as the best 
representation for K of this group. The following numerical values are 
chosen as being representative: 



ap 5 x 10 5 cm^ sec*^ ; a£ .2 x 10 ^ 

* 3 — .99 ; a 4 s .2 x 10 " 3 

The choice of a^ gives a hundred-fold variation of K with height. The 
value of a^ controls the overall magnitude of K: 



For z = 0 , K 5 x 1(P cm^ 6ec”* 



For z=oO , K 5 x 10$ cm^ 



sec 



-1 



The value of a 0 always keeps the value of K positive by limiting the value 



° £ * 2 z y %< 1 



for all elevations and typical lapse rates. The para- 



meter a^ controls the rate of decrease of K with height. 
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3. Solution of Taylor Heat Diffusion Equation . 

The associated boundary conditions of the heat diffusion equation 
(1) are: 



©•= o for z ■ °0 and all t (6) 

&=& 0 (t) for z - 0 (7) 

In the analysis done by Ilaltlner (7) 0 O (t) has been represented 
by two trigonometric terms which is normally a reasonable good approx- 
Lmation to the observed data. Thus the following form is assumed for: 

— 0*5 GJfc -+- c 2 - to fc) 

Here 



' 0 

2.TT 



(8) 



664-0 0 



sec 



-1 



j a .087 which for convenience keeps the 



magnitude ofB 0 5-*l ; c • 0.3 

In order to scale the equation for computational purposes, let 
t » 3600rand z ■ lOOq cr (9) 

Here q is a constant, while <r and T are the new variables, with t 
in seconds, the units of V are hours. For z in cm, cr- will be units 
of q- meters, <-e in meters when q » 1, in 2 meters unit when q - 2, 
etc., with the transformation (9) equations (1) and (4) become: 

(*!<£& -V *£-1 (10) 

h ' v L ^ 3<r ■*«- J 

k . a,[, - <U> 

The form of (8) remains the same except that oO must now be taken 



as 



xTT 

2.4 



. The boundary conditions (6) and (7) remain identical in 



form. It will be worth mentioning here that two-meter level may be 
considered as the lowest level in order to neglect the surface effect. 
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4. Finite Difference Equations. 

In order to obtain a numerical solution of the problem the derivatives 
of 0 with respect to cr in equation (10) and (11) are replaced by appro- 
priate finite difference forms. The problem is then reduced to one of 
solving a system of linear algebraic equations in the values of 0 over a 
grid of points covering the desired range of time and height. We obtain: 



itere the subscript i designates the ith level of the vertical grid at 
height ih, where h is the vertical distance between the points of the grid 
in units of q-meters. Substituting the values of (13) and (14) for i = 1, 2, 



are obtained. Now the partial differential equation (1) has now been reduced 
to a system of ordinary differential equations, which will be integrated 





in the right hand side of (12) the respective values of 




numerically by the Runga - Kutta method. 
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5. Computation. 

3 

For this particular problem q is selected as q «» 10 . This 
means the unit3 of <r are in thousands of meters. In order to 
achieve reasonable accuracy a space interval of ten meters is 
selected in the vertical. This will be sn appropriate time to 
mention that there will be less accuracy in the results below ten 
meters. All the figures below this level are approximated linearly. 
A fairly detailed picture of’ the vertical structure is obtained 
assuming that the upper boundary condition (6) applies at the 
height of 290 meters. Therefore — O . Hence the system 

(iz) to be solved consists of 28 simultaneous ordinary differ- 
ential equations, together with equation (8) for . A numerical 
solution of this system is obtained by the Runga - Kutta method on 
a National Cash Register 102A electronic computer. The choice of 
the appropriate time interval over which the integration is to be 
carried out depends intimately upon the size of the eddy diffusivity. 
As the latter is increased. At must be decreased in order to 
maintain sufficient accuracy. With the value of a^» 5 x 10^ , the 

time interval of the integration period was prohibitively small for 
the type of computer being used, i e., thousands of hours of machine 
time would be required to complete one full cycle. Hence a smaller 

3 

value of a^» 1.92x10 was used, which permitted a time interval of 
1/16 hours. 



12 



% 










6 . 



Results. 



5 3 

The reduction of from 5 x 10 to 1.92 x 10 reduces the 
coefficient of diffusivity which, in turn, reduces the amplitude 
and increases the phase leg with elevation. In Fig. (3) relative 
amplitude and phase lag in minutes are plotted against height in 
meters. It may be observed that the amplitude at ten meters is 
approximately one half of the surface value and the phase lag at 
this level is 78 minutes. These values are not uncommon in winter. 

These values may be compared with those observed at Porton in 
December on clear days. For instance at 17 meters at Porton 
the phase lag is 72 minutes and theory gives 63 minutes. In 
Fig. (3) for comparative purposes, the amplitude reduction and 
phase lag observed at the University of Washington micrometeor- 
logical mast and tower is given in the inset. Notice that in 
the theoretical results amplitude falls off more rapidly at lower 
levels where K is small and the change of phase is smaller than at upper 
levels. The lower values of amplitude and larger phase lags 
compared with observed data once again indicative of the small 

value of K. The value of K in the present analysis at the surface 

1 3 

is K » 1.95 x 10 and at the topmost layer, K * 1.95 x 10 . These 

are very small values particularly for clear days in August on the 

land. 

As given in Table (1), the magnitude of lapse rate (in deg 
per cm) at lower levels, considering above the two meter level, is 
of the order, 10 ; and between ten and forty meters, this order is 

of 10~3. In general the lapse rate continues to decrease with elevation. 



13 



X 








HEIGHT ( METERS) 



* 



I 

I 

f 

i 

r 

3oo r 

^80 

2-6o 

2-4^ 

Z2x> 

2 oo 
l&o 
I 6> O l‘ 
J 4 o 
12.0 : 
loo 
80 - 
Go - 



40 - 
20 



O 



£ 



i 




i 

t 



x- - ~r T-j~ 

fig 3 .. ! . 



4 




i . 

i 

-i 







I 



MEAN R6LATIYE AMPLITUDE AMD PHASE 




i ! , 

• 1 ' 1 • . - 

RELATIVE AMPLITUDE 





14 



■ -r 




~J ’ r T" 



HEaa ft. a. ahd p. l. observed 

, i 

AT UN IV OF WASH MICRCWETEORO- 

* i 

LOGICAL TOWER FOR 23 DAYS OF I 



AUG, 1930 



_i 4- 



RELATIVE AMP - f .. 




l I 



f 








' 













From Che computation, the values at lower levels (about 1U meters) 

-4 -5 

and upper level* are of order 10 and 10 respectively which are 
about one-hundredth of the observed values. This decrease in lapse 
rate may be overcome by selecting a^ in equation (5) hundred times 
large. In equation (12) instead of reducing a^ which decreases the 
magnitude of the coefficient of diffusivity, the factor q may be 
increased to achieve computational stability in the computer. The 

* 

increase in q increases the , finite vertical distance h over which 
the devivatives are evaluated. This implies that the accuracy of 
the finite- difference approximation will decrease. 

As an example of this technique, let q be increased by a factor 
of 10, then K will be equal to 1.95 x 10^, which is a reasonable 
value for the diffusion coefficient. Then the finite vertical dist- 
ance h will be 100 meters. In this case Fig. (3) may be considered 
as plotted with height from zero to 3000 meters. This Interpretation 
of the results may be of value if the temperature changes at higher 
elevations are of primary Interest. 

Fig. (3) also shows the amplitude reduction and phase lag for 
the minimum value of the temperature at various elevations. The 
amplitude reduction for the minimum shows greater reduction than the 
maximum, which apparently corresponds to the lower nocturnal 
diffusivity. There is an appreciable difference between the amplitude 
reduction of the maximum and that of the minimum. For example, at the 
40 meters level, the maximum is about 267. of the surface value of 0 , 
while the minimum is about 197. of the surface minimum. At much lower 
and much higher levels, the difference between amplitude reduction 
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of Che maximum and minimum is le9S. The phase lag increases with 
altitude for both maximum and minimum, but the rate of increase 
decreases with height, because of the increase of the coefficient 
of diffusivity. The phase lag of the minimum is somewhat less 
than that of the maximum which does not seem to be consistent 
with greater amplitude reduction of minima compared to the maxima. 
Perhaps these differences would vary in succeeding maxima and 
minima as a greater time elapsed from the initial conditions. 

Some comparisons between the observed data and theoretical 
results may be made by a proper choice of q which Indicates the 
appropriate diffusion coefficient. Table (2) shows the amplitude 
reduction and phase lag from observation taken at Leafield. The 
agreement here is fair at elevations above 60 meters. The choice 
of q*6 gives the finite difference interval of 60 meters. Table 
(3) shows some observation taken at the University of Washington, 
Seattle. In this case the finite difference interval is 160 
meters, but observations are below this height. However, as a 
first approximation, the amplitude reduction and phase lag are 
interpolated between 0 and 160 meters and tabulated. The results 
fit well, but the accuracy of the interpolation is doubtful. 
Similarly Table (4) compares the theory with observed data at 
Eiffel Tower. Here q is very large giving rise to finite diff- 
erence of h » 600 meters, whereas the heights considered at the 
Eiffel Tower are below 300 meters. Yet the comparison with the 
theory is fair. In Table (5) comparison is made with the results 
taken at Porton. In this case h * 20 meters which is not so great 
as to offset the results of computation. The theory fits well with 
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the observed data. Table (6) compares the amplitude reduction in summer 
and winter at Leafield with theory considering q z 25 and 6 respectively, 
this gives h - 250 and 60 meters respectively. The values of 0 considered 
here are proportional to the ten-meter level. The magnitude of amplitude 
reduction in theory and observation are very close to each other. 



Table. 2. 



Amplitude reduction (A.R.) and phase lag (P.L.) taken for the 
case q- 6 versus observation taken at Leafield, England. 





Theory 




Leafield 




A.R. 


P.L. 


Height 


A.R. 


P.L. 


.73 


33 


25 meters 


.90 


- hours. 


.68 


42 


30 " 


- 


1.20 " 


- 


- 


50 " 


.84 




.52 


1.3 


60 " 


- 


1.5 " 


.43 


1.75 


90 " 


- 


1.66 " 


.41 


1.8 


100 " 


.77 




.37 

» 


2.3 


120 " 


- 




Amplitude reduction and pha 


Table. 3. 
se lag taken 


» 

from case q - 


16 versus 



observation at the University of Washington, Seattle. 





Theory 




Observation 


A.R. 


P.L. 


Height 


A.R. 


P.L. 


.91 


12 min 


15 meters 


.91 


16 min 


.83 


18 " 


30 " 


.86 


18 " 


.79 


27 " 


45 " 


.83 


23 " 


.73 


33 " 


60 " 


- 


_ it 
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Table .4. 



Amplitude reduction from Eiffel tover versus theoretical values 
for q - 60. 



Theory 




Observation 




A.R. 


Height (meters ) 


Sunnier 


W inter 


.79 


195 


. 88 


.83 


.68 


300 


.75 


.60 


Table. 5. 

Amplitude reduction and phase lag taken for 


the case q - 2 


versus 


obsertion taken at Porton, England. 






Theory 




Porton 




A.R. P.L. 


Height 


A.R. 


P.L. 


.75 .57 


7 meters 


.73 


.95 hours 


.58 1.05 


17 " 


.65 


1.2 " 




Table. 6. 






Amplitude reduction from Leafield (England) 


versus theoretical values. 


Height in meters 


10 25 50 75 


100 


Relative amplitude. 








(a) June. 


1.00 0.90 0.84 0.79 


0.77 


(b) December. 


1.00 0.74 0 


.54 0.44 


0.40 


Theory: q - 25 


1.00 0.90 0.84 0.78 


0.64 


q - 6 


1.00 0.64 0 


.58 0.47 


0.42 
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7 . 



Conclus ion* 



A numerical solution has been presented for the diurnal temper- 
ature variation with a coefficient of eddy diffusivlty which is a 
function of height and lapse rate. By suitable selection of several 
parameters, reasonably good agreement has been obtained between 
theory and observation. 

Improvement in the results may possibly be achieved by improving 
the functional form of diffusivlty coefficient. It would be particul- 
arly desirable to make a detailed study of the relationship between 
eddy diffusivlty and such parameters as surface roughness, wind speed, 
in addition to stability. If a suitable relationship could be obtained 
for K in terms of these parameters then the solution for the diurnal 
temperature wave would have a broader application. 
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